{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "wicked-shoulder",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "  0%|          | 0/360 [00:00<?, ?it/s]\u001b[A\u001b[A\u001b[A\u001b[A"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "burning.......\n",
      "sampling: iterations 129 of 10000, elapsed 0:00:06, remaining 0:07:03\n",
      "sampling: iterations 337 of 10000, elapsed 0:00:11, remaining 0:05:15\n",
      "sampling: iterations 599 of 10000, elapsed 0:00:16, remaining 0:04:18\n",
      "sampling: iterations 785 of 10000, elapsed 0:00:22, remaining 0:04:13\n",
      "sampling: iterations 972 of 10000, elapsed 0:00:27, remaining 0:04:07\n",
      "sampling: iterations 1159 of 10000, elapsed 0:00:32, remaining 0:04:01\n",
      "sampling: iterations 1459 of 10000, elapsed 0:00:40, remaining 0:03:55\n",
      "sampling: iterations 1719 of 10000, elapsed 0:00:47, remaining 0:03:44\n",
      "sampling: iterations 2020 of 10000, elapsed 0:00:55, remaining 0:03:36\n",
      "sampling: iterations 2281 of 10000, elapsed 0:01:01, remaining 0:03:27\n",
      "sampling: iterations 2585 of 10000, elapsed 0:01:09, remaining 0:03:19\n",
      "sampling: iterations 2848 of 10000, elapsed 0:01:16, remaining 0:03:10\n",
      "sampling: iterations 3152 of 10000, elapsed 0:01:24, remaining 0:03:03\n",
      "sampling: iterations 3415 of 10000, elapsed 0:01:30, remaining 0:02:54\n",
      "sampling: iterations 3719 of 10000, elapsed 0:01:39, remaining 0:02:47\n",
      "sampling: iterations 3982 of 10000, elapsed 0:01:45, remaining 0:02:39\n",
      "sampling: iterations 4286 of 10000, elapsed 0:01:54, remaining 0:02:32\n",
      "sampling: iterations 4475 of 10000, elapsed 0:01:59, remaining 0:02:27\n",
      "sampling: iterations 4664 of 10000, elapsed 0:02:04, remaining 0:02:22\n",
      "sampling: iterations 4927 of 10000, elapsed 0:02:10, remaining 0:02:13\n"
     ]
    }
   ],
   "source": [
    "from data_generation import Long_Data_gener\n",
    "import numpy as np\n",
    "from scipy import special\n",
    "from nonparameter_classfication import NonparaClassfication\n",
    "import parameter_classfication\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "\n",
    "import pyjags\n",
    "def dina_parameter_estimate(N,Q,I,K,simu_resp,chains=2,cores=4,burning=5000,iteration=5000,return_trace=False,return_full_vars=False):\n",
    "    # hodina_pyjags\n",
    "    jags_ho_dina_data = {\n",
    "        \"Y\": simu_resp,\n",
    "        \"Q\": Q,\n",
    "        \"N\": N,\n",
    "        \"K\": K,\n",
    "        \"I\": I\n",
    "\n",
    "    }\n",
    "    if return_full_vars == True:\n",
    "        jags_model_return_vars = [\"p\", \"alpha_p\", \"alpha\", \"theta\", \"beta\", \"lambda_\", \"s\", \"g\"]\n",
    "    else :\n",
    "        jags_model_return_vars = [\"alpha\"]\n",
    "\n",
    "    jags_hodina_model_str =\\\n",
    "    \"\"\"\n",
    "    model{\n",
    "        for (n in 1:N) {\n",
    "            for (i in 1:I){\n",
    "                for (k in 1:K){\n",
    "                    w[n,i,k] <-(alpha[n,k]>=Q[i,k])\n",
    "                }\n",
    "                eta[n,i] <- prod(w[n,i,1:K])\n",
    "                p[n,i] <- eta[n,i]*(1-s[i]-g[i]) + g[i]\n",
    "                Y[n,i] ~ dbern(p[n,i])\n",
    "            }\n",
    "        }\n",
    "\n",
    "        for (n in 1:N){\n",
    "            for (k in 1:K){\n",
    "                logit(alpha_p[n,k]) <- lambda_[k]*theta[n]+beta[k]\n",
    "                alpha[n,k] ~ dbern(alpha_p[n,k]) \n",
    "            }\n",
    "            theta[n] ~ dnorm(0,1)\n",
    "        }\n",
    "\n",
    "        for (k in 1:K){\n",
    "            beta[k] ~ dnorm(0,0.25)\n",
    "            lambda_[k] ~ dnorm(0,0.25) T(0,)\n",
    "        } \n",
    "        for (i in 1:I){\n",
    "            s[i] ~ dbeta(1,1)\n",
    "            g[i] ~ dbeta(1,1) T(,1-s[i])\n",
    "        }\n",
    "    }\"\"\"\n",
    "    jags_ho_dina_model = pyjags.Model(code=jags_hodina_model_str,\n",
    "                                        data=jags_ho_dina_data,\n",
    "                                        chains=chains,\n",
    "                                        threads=cores,\n",
    "                                        generate_data=True,\n",
    "                                        chains_per_thread=1,\n",
    "                                        adapt=False,\n",
    "                                      progress_bar=True\n",
    "                                        )\n",
    "    print(\"burning.......\")\n",
    "    jags_tarce_prior = jags_ho_dina_model.sample(iterations=burning, vars=[], thin=1)\n",
    "    print(\"sampling.......\")\n",
    "    jags_tarce_posterior = jags_ho_dina_model.sample(iterations=iteration, vars=jags_model_return_vars, thin=1)\n",
    "    x = (jags_tarce_posterior[\"alpha\"].mean(axis=-1).mean(axis=-1) > 0.5)\n",
    "    if return_trace==True:\n",
    "        return x,jags_tarce_posterior\n",
    "    return x\n",
    "\n",
    "\n",
    "def long_dina_parameter_estimate(T,N,Q,I,K,simu_resp,burning=5000,iteration=5000,chains=2,cores=4\n",
    "                                 ,return_trace=False,return_full_vars=False):\n",
    "    jags_data = {\n",
    "        \"Y\": simu_resp,\n",
    "        \"Q\": Q,\n",
    "        \"N\": N,\n",
    "        \"I\": I,\n",
    "        \"K\": K,\n",
    "        \"T\": T\n",
    "    }\n",
    "    if return_full_vars==True:\n",
    "        jags_return_vars = [\"mu_theta\", \"pr_theta\", \"theta\", \"att_possi\",\n",
    "                        \"att\", \"theta_beta\", \"theta_lambda\", \"lambda_0\", \"lambda_k\", \"p\"]\n",
    "    else :\n",
    "        jags_return_vars = [\"att\",\"att_possi\"]\n",
    "\n",
    "    long_dina_jags_model = ling_dina_jags_model = \\\n",
    "        \"\"\"\n",
    "        model{\n",
    "            for (t in 1:T){\n",
    "                for (n in 1:N){\n",
    "                    for (i in 1:I){\n",
    "                        for (k in 1:K){w[t,n,i,k] <- (att[t,n,k]>=Q[t,i,k])}\n",
    "                        eta[t,n,i] <- prod(w[t,n,i,1:K])\n",
    "                        logit(p[t,n,i]) <- lambda_0[t,i]+lambda_k[t,i]*eta[t,n,i]\n",
    "                        Y[t,n,i] ~ dbern(p[t,n,i])\n",
    "                    }\n",
    "                }\n",
    "    \n",
    "    \n",
    "                for (n in 1:N){\n",
    "                    for (k in 1:K){\n",
    "                        logit(att_possi[t,n,k]) <- theta_beta[k] + theta_lambda[k]*theta[t,n]\n",
    "                        att[t,n,k] ~ dbern(att_possi[t,n,k])\n",
    "                    }\n",
    "                }\n",
    "            }\n",
    "    \n",
    "            for (n in 1:N){\n",
    "                theta[1:T,n] ~ dmnorm(mu_theta[1:T], pr_theta[1:T,1:T])\n",
    "            }\n",
    "            for (k in 1:K){\n",
    "                theta_beta[k]~ dnorm(0,0.25)\n",
    "                theta_lambda[k] ~ dnorm(0,0.25) T (1,)\n",
    "            }\n",
    "    \n",
    "            for (t in 1:T){\n",
    "                for (i in 1:I){\n",
    "                    lambda_0[t,i] ~ dnorm(-1.096, 0.25)\n",
    "                    lambda_k[t,i] ~ dnorm(0, 0.25) T(0, )\n",
    "                }\n",
    "            }\n",
    "    \n",
    "            mu_theta[1] <- 0\n",
    "            for (t in 2:T){\n",
    "                mu_theta[t] ~ dnorm(0,0.5)\n",
    "            }\n",
    "    \n",
    "            L_theta[1,1] <- 1\n",
    "            for (tt in 2:T){\n",
    "                L_theta[tt,tt] ~ dgamma(1,1)\n",
    "                for (ttt in 1:(tt-1)){\n",
    "                    L_theta[tt,ttt] ~ dnorm(0,1)\n",
    "                    L_theta[ttt,tt] <- 0\n",
    "                }\n",
    "            }\n",
    "            sigma_theta = L_theta %*% t(L_theta)\n",
    "            pr_theta[1:T, 1:T] <- inverse(sigma_theta[1:T, 1:T])\n",
    "    \n",
    "    \n",
    "        }\n",
    "        \"\"\"\n",
    "    jags_long_dina_model = pyjags.Model(code=long_dina_jags_model,\n",
    "                                        data=jags_data,\n",
    "                                        chains=chains,\n",
    "                                        threads=cores,\n",
    "                                        generate_data=True,\n",
    "                                        chains_per_thread=1,\n",
    "                                        adapt=False,\n",
    "                                        progress_bar=True\n",
    "                                        )\n",
    "    print(\"burning.......\")\n",
    "    jags_tarce_prior = jags_long_dina_model.sample(iterations=burning, vars=[], thin=1)\n",
    "    print(\"sampling.......\")\n",
    "    jags_tarce_posterior = jags_long_dina_model.sample(iterations=iteration, vars=jags_return_vars, thin=1)\n",
    "\n",
    "    x = (jags_tarce_posterior[\"att\"].mean(axis=-1).mean(axis=-1) > 0.5)\n",
    "    if return_trace==True:\n",
    "        return x,jags_tarce_prior\n",
    "    return x\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def get_indicator(att_pre, att_real):\n",
    "    x = att_pre == att_real\n",
    "    ACCR = np.nanmean(x)\n",
    "    y = np.product(x, axis=-1)\n",
    "    PCCR = np.nanmean(np.product(x, axis=-1))\n",
    "    Long_PCCR = np.nanmean(np.product(np.product(x, axis=-1), axis=0))\n",
    "    return ACCR, PCCR, Long_PCCR\n",
    "\n",
    "\n",
    "def run(t, k, i, N, s, g):\n",
    "    theta_beta = np.ones((k)) * 1.5\n",
    "    theta_lambda = np.ones((k)) * (-1) + np.array([i * 0.5 for i in range(k)])\n",
    "    lambda_0 = np.repeat(special.logit(g), i * t).reshape(t, i)\n",
    "    lambda_k = np.repeat(special.logit(1 - s) - special.logit(g), i * t).reshape(t, i)\n",
    "    # 实例化对象\n",
    "    gener = Long_Data_gener()\n",
    "    classf = NonparaClassfication()\n",
    "    # 生成被试高阶能力\n",
    "    theta_cov = np.array([[1, 1.25, 1.76, 2.74, 4.26],\n",
    "                         [1.25, 1.56, 1.76, 2.19, 2.73],\n",
    "                         [1.76, 1.76, 2.43, 3.8, 4.26],\n",
    "                         [2.74, 2.19, 2.73, 3.8, 4.26],\n",
    "                         [4.26, 2.73, 3.41, 4.26, 5.9]])\n",
    "    if t == 5:\n",
    "        theta = gener.m_theta(t, mean=0, growth_rate=0.5, n=N, input_cor_mat=theta_cov)\n",
    "    elif t==3:\n",
    "        theta = gener.m_theta(t, mean=0, growth_rate=0.5, n=N, input_cor_mat=np.array(theta_cov[:3,:3]))\n",
    "    elif t==2:\n",
    "        theta = gener.m_theta(t, mean=0, growth_rate=0.5, n=N, input_cor_mat=np.array(theta_cov[:2,:2]))\n",
    "\n",
    "    # 生成被试属性\n",
    "    att, z = gener.attribute(theta, theta_beta, theta_lambda, return_possi=True)\n",
    "    # 生成题目Q矩阵\n",
    "    Q = gener.Q_generate_long(t, i, k)\n",
    "    # 生成模拟作答\n",
    "    resp, resp_possi = gener.simu_resp(t, att, Q, lambda_k=lambda_k, lambda_0=lambda_0, return_possi=True)\n",
    "\n",
    "    \"\"\"被试数据生成到此结束\"\"\"\n",
    "    # 获取被试所有可能的属性掌握模式\n",
    "    all_pattern = gener.all_pattern(k)\n",
    "    # 计算所有可能的属性掌握模式对应的理想作答\n",
    "    class_ideal_resp = gener.long_dina_ideal_resp(t, Q, all_pattern)\n",
    "    # 参数方法\n",
    "    para_pre_att = dina_parameter_estimate(N,Q[0],I=i,K=k,simu_resp=resp[0],return_trace=False)\n",
    "    # para_pre_att = att[0]\n",
    "    # para_pre_att = att[0]\n",
    "    # 使用非参方法对被试进行分类\n",
    "    class_result_dic = {}\n",
    "    class_result_dic[\"class_result1\"] = classf.long_standard_distance(t, all_pattern, simu_resp=resp,\n",
    "                                                                      class_ideal_response=class_ideal_resp,\n",
    "                                                                      distance=\"hamming\")\n",
    "    class_result_dic[\"class_result2\"], l = classf.long_standard_dis_constrain_1(t, all_pattern, simu_resp=resp,\n",
    "                                                                                class_ideal_response=class_ideal_resp,\n",
    "                                                                                distance=\"hamming\", return_judge=True)\n",
    "    class_result_dic[\"class_result3\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "                                                                             class_ideal_response=class_ideal_resp,\n",
    "                                                                             distance=\"hamming\")\n",
    "    class_result_dic[\"class_result4\"], l = classf.long_standard_dis_constrain_1(t, all_pattern, simu_resp=resp,\n",
    "                                                                                class_ideal_response=class_ideal_resp,\n",
    "                                                                                distance=\"hamming\", return_judge=True,\n",
    "                                                                                input_first_t_att=att[0, :, :])\n",
    "    class_result_dic[\"class_result5\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "                                                                             class_ideal_response=class_ideal_resp,\n",
    "                                                                             distance=\"hamming\",\n",
    "                                                                             input_first_t_att=att[0, :, :])\n",
    "    class_result_dic[\"class_result6\"] = classf.long_standard_distance(t, all_pattern, simu_resp=resp,\n",
    "                                                                      class_ideal_response=class_ideal_resp,\n",
    "                                                                      distance=\"hamming\", weight=True)\n",
    "    class_result_dic[\"class_result7\"], l = classf.long_standard_dis_constrain_1(t, all_pattern, simu_resp=resp,\n",
    "                                                                                class_ideal_response=class_ideal_resp,\n",
    "                                                                                distance=\"hamming\", return_judge=True,\n",
    "                                                                                input_first_t_att=para_pre_att)\n",
    "    class_result_dic[\"class_result8\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "                                                                             class_ideal_response=class_ideal_resp,\n",
    "                                                                             distance=\"hamming\",\n",
    "                                                                             input_first_t_att=para_pre_att)\n",
    "    class_result_dic[\"class_result9\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "                                                                             class_ideal_response=class_ideal_resp,\n",
    "                                                                             distance=\"weight_hamming\",\n",
    "                                                                             input_first_t_att=att[0, :, :])\n",
    "    class_result_dic[\"class_result10\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "                                                                              class_ideal_response=class_ideal_resp,\n",
    "                                                                              distance=\"weight_hamming\",\n",
    "                                                                              input_first_t_att=para_pre_att)\n",
    "    # class_result_dic[\"class_result11\"] = classf.long_standard_dis_constrain_2(t, all_pattern, simu_resp=resp,\n",
    "    #                                                                           class_ideal_response=class_ideal_resp,\n",
    "    #                                                                           distance=\"weight_hamming\",\n",
    "    #                                                                           input_first_t_att=para_pre_att)\n",
    "    class_result_dic[\"class_result11\"] = long_dina_parameter_estimate(T=t,N=N,Q=Q,I=i,K=k,simu_resp=resp)\n",
    "\n",
    "    class_mathod_name_list = [\"标准\", \"限制1\", '限制2',\n",
    "                              '第一次完全正确限制1', '第一次完全正确限制1',\n",
    "                              \"加权距离\", \"第一次为参数估计限制1\", \"第一次为参数估计限制2\",\n",
    "                              \"第一次完全正确加权距离限制2\", \"第一次为参数估计加权距离限制2\", \"Long_dina_参数\"]\n",
    "    rename_dict = {}\n",
    "    for i in range(len(class_mathod_name_list)):\n",
    "        att_pre = class_result_dic[\"class_result{}\".format(i + 1)]\n",
    "        A, P, L = get_indicator(att_pre, att)\n",
    "        rename_dict[\"ACCR_\" + class_mathod_name_list[i]] = A\n",
    "        rename_dict[\"PCCR_\" + class_mathod_name_list[i]] = P\n",
    "        rename_dict[\"Long_PCCR_\" + class_mathod_name_list[i]] = L\n",
    "\n",
    "    return rename_dict\n",
    "\n",
    "\n",
    "R = 15\n",
    "I_list = [50, 30]\n",
    "T_list = [5,2,3]\n",
    "K_list = [3, 5]\n",
    "N_list = [300, 100]\n",
    "\n",
    "pbar = tqdm(total=len(I_list) * len(T_list) * len(K_list) * len(N_list) * R)\n",
    "for t in T_list:\n",
    "    for k in K_list:\n",
    "        for N in N_list:\n",
    "            for i in I_list:\n",
    "                result_list = []\n",
    "                for r in range(R):\n",
    "                    rename_dict = run(t=t,\n",
    "                                      k=k,\n",
    "                                      i=i,\n",
    "                                      N=N,\n",
    "                                      s=0.1,\n",
    "                                      g=0.1)\n",
    "                    result_list.append(rename_dict)\n",
    "                    pbar.update(1)\n",
    "\n",
    "                df = pd.DataFrame(result_list)\n",
    "                df.to_csv(\"./result2/result_t{}_k{}_n{}_i{}.csv\".format(t, k, N, i), index=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "restricted-navigation",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "incorrect-playback",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "patent-progress",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
